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Summary of Technical Sessions and Discussions 

A MECA workshop, “Dust on Mars D,” was sponsored by NASA through the Lunar and Planetary 
Institute and hosted by Arizona State University on February 24-25, 1986. Following the recommendations 
of the participants in the previous MECA “dust” workshop (held at ASU on February 4-5, 1985), the 
goal of this workshop was to discuss the following questions: 

(1) How many components of dust are there on Mars, and what are their properties? 

(2) How is dust ejected from the surfeice into the atmosphere? 

(3) How does the global atmospheric circulation affect the redistribution of dust? 

(4) Are there sources and sinks of dust? If so, how do they vary with time? 

Forty-three people attended the meeting and engaged in the discussion sections; of these, twenty-four 
presented summaries of their dust-related research. 

The first workshop session considered the physical and chemical properties and the distribution (both 
temporal and spatial) of dust and condensates. Several “observational” presentations focused on detailed 
analyses of Viking data. Pressure data recorded over four martian years at Viking Lander 1 and two 
years at Viking Lander 2 have been studied by J. Tillman with particular attention to the major transient 
events related to two global dust storms in 1977 (Year 1) and one in 1982 (Year 4). The 1982 event 
was similar in season to the first 1977 storm (L, ~200°) but was more intense than this storm, appearing 
instead to be similar to the second 1977 storm. The 1982 storm was noted as having a unique pressure 
signature, suggesting a second surge of intense activity a week after the initial change in pressures was 
noted. T. Martin presented a dust opacity history based on a global data analysis using the 7, 9, and 
15 /im data collected by the Orbiters’ Infrared Thermal Mappers (IRTM). The resultant maps of opacity 
depict the spatial and temporal variations of atmospheric dust loading throughout the Viking missions. 
Detailed maps during the 1977 dust storm season suggest large local variability of opacity in many areas. 

Evidence for HjO ice in the Mars atmosphere has been inferred from a variety of optical observations 
using the Landers and Orbiters. Limb observations show that a detached haze fluctuates in height and 
opiacity with time and location. R. Kahn modelled the properties of these hazes during the non-dust storm 
periods with consideration of mean optical depth, condensation level, mean particle size, and eddy diffusion. 
Correlation of the cloud height with MAWD column water vapor abundance implies that water controls 
the location of the cloud base. The mean particle sizes increase as condensation heights decrease in 
such a way that the average fall time is always about 1/4 day. F. Jaquin reported on analysis of the 
limb observations made by the Orbiter cameras over three Mars years. Profiles of reflectance vs. elevation 
above the surface show the buildup of dust with height as a function of time and latitude. During the 
1977 dust storm season, continuous haze up to 50 km altitude is seen in the northern hemisphere with 
a weak detached haze above that. The detached haze (probably condensates) is better defined in the 
southern hemisphere. 

The second set of speakers reported on “laboratory emd theoretical” results. J. Gooding has been 
performing experiments and modelling the heterogeneous nucleation of condensates on mineral grains. 
Favorable conditions for condensation appear to be related to crystallographic compatibility (with respect 
to structure and number of “active” sites) and the nature of chemical bonds. Details for mineral classes 
and for given minerals, especiaUy clays, show the complexities of the nucleation process as a function 
of condensate type. Furthermore, as nucleation occurs the initial condensates will change conditions 
associated with subsequent nucleation. D. Colburn discussed the influence of dust on water condensation 
using Lander optical depth measurements. Both seasonal and diurnal effects appear to be factors in the 
amount of water condensation. Also, most condensation appears to occur aloft as opposed to near-surface, 
and is strongly related to dust abundance. Simulating results from the labeled-release biology experiment 
aboard the Landers, A. Banin inferred properties of martian dust. Smectite minerals (montmorillonite and 
nontronite) that were ion-modified were found to best simulate the kinetic results of the biology experiment. 
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The ion-modification served to exchange calcium and sodium ions with iron. This modification is consistent 
with Lander chemical analyses. These experiments suggest that: (1) the pH of martian soil is less than 
4.5; (2) CaCOs probably cannot be in soil at concentrations larger than ~1.0%; (3) montmorillonite is 
apparently essential for causing rapid decomposition on clay surfaces. According to results presented by 
R. Morris, the spectra of martian bright regions suggest that "dust” produced in the dry valleys of Antarctica 
may be a reasonable martian analogue. Weathering on Mars, assuming analogous processes to those 
in Antarctica, results in more highly oxidized material, but in general appears to be similar to that for 
ferric products that were studied. Further, the presence of hematite as a component of martian dust 
is strongly suggested. 

The second workshop session addressed topics of atmospheric dust transport and redistribution. L. 
Martin led off the session discussing equatorial dust clouds imaged by Viking that differed from other 
dust clouds because of their persistent association with a small canyon. Downslope katabatic-type winds 
may be responsible for the clouds’ oryn; in some cases they appeared to flow into the canyon floor. 
The potential of the telescopic record for revealing interannual variation in the CO 2 cycle was the subject 
of the presentation $ven by P. James. The signature of this variability appears as (Merences in the retreat 
rate of the south polar cap (readily observed during favorable oppositions). Comparison of data for the 
1956, 1971, 1973, and 1977 dust storm years shows some variability with the ca^} retreating slower than 
the median in 1971 and 1977 and faster than the median in 1956 and 1973. The extent to which this 
variability would affect pressure variations is evidently small, however, given the high degree of repeatability 
in the seasonal pressure data from the Lander sites. D. Paige discussed his analysis of IRTM data at 
the north and south poles during the spring season. By constraining a radiative equilibrium model he 
developed to reproduce the brightness temperatures in the various IRTM channels, he concluded that 
dust in both polar regions is similar in composition to that observed in the tropics (the so called “Toon- 
Pollack” dust). He also found that unless the dust was confined to the lowest 5 km, it was not possible 
to reproduce the observed 15 micron brightness temperatures. This problem was even more difficult in 
the north where additional cooling, perhaps due to water ice, is needed. 

J. Barnes focused attention on the transport implications of the warming of the north polar atmosphere 
during the second global dust storm of 1977. Based on models that follow the actual motion of air parcels, 
the observed warming must be accompanied by sinking motion; to conserve mass, poleward transport 
from lower latitudes is required. Barnes also discussed other modelling that implies that vertically-propagating 
forced stationary waves may play an important role in the heat budget of the polar atmosphere during 
global dust storms. In the next paper, R. Haberle presented a mechanism for interannual variability of 
global dust storms. His nearly inviscid zonally-symmetric circulation model suggests the following scenario: 
First, global dust storms transport dust into the northern hemisphere. This dust is then available for raising 
by baroclinic waves in the following years. During the post-dust-storm years the cross-equatorial Hadley 
circulation is therefore weakened due to heating by dust suspended in its descending branch. When the 
supply of dust is exhausted, the full strength of the Hadley circulation, and the likelihood of global dust 
storms, is restored. R. Zurek then spoke about how tidal forcing of the mean meridional circulation may 
break up the otherwise smooth structure of the cross-equatorial Hadley cell. During relatively dusty periods, 
tides produce significant accelerations of the zonal and meridional wind fields, and these accelerations 
have considerable vertical structure. When included in Haberle’s 2D model, several vertically-stacked Hadley- 
type cells were produced. Zurek suggested that by changing the structure of the Hadley circulation, the 
tides could effectively limit the flux of dust into the northern hemisphere during a global dust storm, 
explaining the decay patterns seen at the Lander sites. Completing the “theoretical” presentations of this 
session, J. Pollack gave a progress report on the Mars general circulation model (GCM) and presented 
some preliminary results. Work on the GCM during the past several years has focused on including dust 
as a radiativcly active constituent, generalizing the model to an arbitrary number of layers, and improving 
the numerical methods for calculating infrared fluxes. All of these changes have been implemented and 
tested, and the current effort has shifted to determining the model’s performance in a variety of configurations 
and examining how dust affects the behavior of large-scale atmospheric eddies. 


During the final session of the workshop, the participants turned their attention to the study of dust 
sedimentation and erosion on Mars. The research that was discussed is being conducted in three ^neral 
areas: photointerpretation of sedimentary and erosional features, surface properties determined from thermal 
inertia, albedo, and radar data, and laboratory and theoretical work on the processes involved in wind 
transport of dust. R. Greeley tied together the various laboratory and theoretical studies applicable to 
the martian dust cycle. In particular, the continuing question of the role of volatiles in any aggregation 
of dust particles and in ejection of dust from the surface was considered. R. Arvidson presented results 
of a study of the morphology and statigraphy of south polar pitted and layered terrains that indicates 
deposition in two very different phases extending over about half of martian history. The radically different 
layered and unlayered pitted deposits suggest substantial differences in sedimentary regimes. At another 
size scale, H. Moore reported that the morphology of small landslides at the VLl site is diagnostic of 
surface layers of. very low cohesiveness that are probably recent aeolian deposits. Slope failure of ^gch 
layers may be an important mechanism for producing a surface that is not in equilibrium with maraan 
wind regimes, thereby enhancing the possibility of aeolian erosion. 

Several specific regions of the planet were discussed in detail. The physical properties of aeolian features 
in the Elysium-Amazonis area were examined by J. Zimbelman. Analysis of the thermal inertia of wind 
streaks in the region suggests a long-lived movement of fine to medium sand-sized materials in the observed 
depositional dark streaks and possible dust deposits of significant thickness accumulating in bright streaks. 
B. Jakosky reported on a study relating the characteristics of the Lander sites to the global data on 
thermal inertia, albedo, and color; the sites apparently do not match the general planet-wide trends in 
these data sets. In particular, the VLl site is intermediate between bright areas assumed to accumulate 
dust and those dark areas that are thought to be swept relatively clear. P. Christensen discussed the 
albedo variation of several areas on the planet using the IRTM albedo channel. He found that northern 
hemisphere dark regions showed a long period of progressive dust removal after global dust storms, while 
southern hemisphere dark regions show a rapid return to pre-storm albedos. S. Lee investigated regional 
sources and sinks using IRTM albedo data and Orbiter images. The patterns of temporal changes indicate 
that Syrtis Major, through much of the year, acts as a source region for dust redeposited in neighboring 
Arabia. In the Solis Planum region, removal of dust is restricted to periods of major dust-storm activity. 

The ability of radar to detect seasonal variability in surface dust cover was examined by L. Roth. A 
theoretical study of radar reflectivity as a function of dust layer thickness concludes that time-variable 
reflectivity data would be consistent with significant changes in dust layer thickness only under very restricted 
conditions of layer geometry and continuity. The session was concluded with T Thompson’s review of 
future opportunities for radar observations of Mars. The 1986 and 1988 oppositions will greatly extend 
the southern hemisphere coverage available from previous years, while upgrades in the Goldstone radar 
should significantly improve resolution and allow both X- and S-band observations. 

The final afternoon of the workshop consisted of an open discussion of the major questions posed 
during the preceding two days. Much of the discussion centered around the possibilities for interdisciplinary 
cooperation. For example, observations of regional dust sources and sinks, and tracking the location, 
timing, and growth of local and great dust storms, should provide useful input and constraints for atmospheric 
circulation models. The mineralogical properties of dust may greatly influence the formation of condensates; 
hence, knowledge of such properties may play major roles in modelling the production of atmospheric 
hazes and polar deposits. It was agreed that defining the frequency and characteristics of local and “global” 
dust storms, the distribution (temporal and spatial) of dust on the surface, and the physical and chemical 
properties of dust along with inclusion of dust related effects in atmospheric circulation models, are all 
areas in which continued research and communication among the various interested groups would be 
fruitful for the entire community. The research discussed at the workshop is clearly moving in the direction 
of posing questions for the Mars Observer mission, while showing that Viking Orbiter and Earth-based 
data provide an excellent basis for formulation of Mars Observer experiment strategies. 
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ON THE RATE OF FORMATION OF SEDIMENTARY DEBRIS ON MARS; Raymond 
E. Arvidson, McDonnell for the Space Sciences, Department of Earth and 
Planetary Sciences, Washington University, St. Louis, MO 63130 

We now have enough information to place rather meaningful constraints on 
the current rates of aeolian redistribution of fine grained debris on Mars. 
For example, from the three years of Viking Lander 1 observations, typical 
values of sediment redistribution are micrometers per year, although 
centimeters of loose, disturbed material were removed during a brief interval 
in the third year (1). Rates of erosion of rocks were too small to be 
observed, either by tracking changes in morphologies or by tracking changes 
in rock photometric properties. Consideration of the large number of 
pristine-looking, small bowl-shaped craters at the Lander 1 site suggests a 
rate of rock breakdown and removal of only meters over the lifetime of the 
surface. Thus, averaged over the lifetime of the Chryse Plains, rock 
breakdown and removal has been meters per billions of years, orders of 
magnitude lower than the micrometers per year for soils (2). Most of the 
equatorial regions of Mars likewise preserves ancient surfaces, with craters 
even at small sizes seemingly in production. Thus, rock breakdown and 
removal over the equatorial terrains has been quite small for much of 
geological time. Even the fretted terrains have probably been inactive for a 
long while, considering that the crater abundances in much of fretted terrain 
are second in abundance only to the cratered terrains (3). The younger 
fretted areas seem to be embayed by younger deposits (4) or to be in areas of 
relatively recent tectonic activity, such as the chaotic terrains. By areal 
extent, most of the equatorial terrains of Mars have been subjected to very 
low erosion rates, significantly less than the debris redistribution rates 
witnessed by Viking Lander 1. Thus, as noted by (2), differential aeolian 
erosion on Mars is a major geological process, with debris deposits 
accumulating and eroding to depths of hundreds of meters over geological 
time, while rock breakdown has been occurring in most regions at vanishingly 
small rates. Given the low rates of production of new debris, one is forced 
to conclude that Mars has had a decidedly non-linear history of debris 
production. In particular, most sedimentary debris must have been produced 
relatively early, perhaps in the first billion years of geological time. 
Impact cratering, production of volcanic ash, and chemical corrosion may all 
have been important debris-forming processes. Given that the mineralogy of 
martian debris has apparently not come into equilibrium with the present 
ambient conditions (5), we may have some chance of deciphering the relative 
importance of various debris forming processes. For example, if analyses of 
Mars Observer Imaging Spectrometer data show a dominance of palagonitic 
materials, then early volcanism would be indicated. 

REFERENCES 
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CONSTRAINING MARS-DUST MINERALOGY ON THE BASIS OF VIKING BIOLOGY 
SIMULATIONS AND MARS SPECTRAL REFLECTANCE. 

A. Banin, Hebrew University, P.O. Box 12, Rehovot, Israel 

Among the prime, unresolved and often debated issues related to Mars is 
the question of the composition and properties of its fine, dusty surface 
materials. Unfortunately, no direct and definitive mineralogical analysis has 
been done on the soil and dust of Mars as yet. Various indirect approaches 
were used, however, to arrive at Mars dust mineralogy. These includ^ direct 
chemical analyses, chemical modeling, simulations of the chemical interaction 
with water and solutes, and a multitude of spectroscopic observations from 
Earth, orbiters and landers. In the following we will briefly review these, 
specifically emphasizing the chemical simulations and present our current 
thinking on the mineral composition of the Martian dust. 

Chemical composition of the Martian fine surface material; Quite 
probably, the most directly related available evidence regarding soil and dust 
mineralogy, is the inorganic chemical analysis (ICA) obtained by X-ray 
fluorescence analyzers aboard the Viking landers. The analyses have shown 
elemental abundances, (given as oxides), of 40-45% Si02, 7-7.5% AI2O3, 17-19% 
Fe20 3, 5-7% MgO, 5-6% CaO, 0.4-0.7% Ti02, traces of K, 6-9% SCh and 0. 3-0.9% 
Cl (1). The chemical composition of the soil in the two Viking landing sites, 
and of different samples taken from each site, was found to be strikingly 
similar. On the basis of the chemical data, several combinations of minerals 
have been suggested initially (2,3,4). These characteristically contained 
about 60-80% clay minerals of the smectite group, mixed with various soluble 
salts such as kieserite (MgSO^ ) and halite (NaCl). Suggested accessory 
minerals were: iron oxides(Fe 0 ), calcite (CaC03), and quartz (Si02). 

Simulation of the Viking^ Biology experiments: The Viking Biology 
experiments turned out to be (and still are) an important and unique source of 
chemical information on Martian dust and soil because of the direct "wet 
chemistry" that was conducted on the soil and the information obtained on the 
interaction of the Mars minerals with water, salts and organics. Although no 
clear evidence for life was found, the Viking Biology experiments showed that 
the minerals in the Mars soil were chemically reactive and capable of 
decomposing various organic acids, catalyzing photochemical fixation of 
CO2/CO, and releasing oxygen upon wetting (5,6). 

The simulation of the Viking Biology experiments, conducted between 1977 
and 1979 gave somewhat diverse results. Levin and his colleagues reported 
that the Mars soil analogs provided by the ICA Team did not simulate the 
Labeled Release (LR) Viking experiment (7). On the other hand, Banin and his 
colleagues, using ion-modified clays, succeeded in the kinetic simulations of 
this experiment (8,9,10). The clay modification by Banin et al. (11,12) 
involved the exchange of calcium and sodium ions, the primary adsorbed ions of 
clays on Earth, with iron, a ubiquitous and abundant element in the Martian 
soil. In addition to affecting the chemical and catalytic properties of the 
clay, this exchange resulted in a pronounced color change from gray to 
reddish-orange, making the clay very similar in color and spectral 
characteristics to the Martian soil (13,14). 

Our detailed experimental studies showed (10,13,14,15) that among the 
many mixtures of clays and various salts whose combinations give elemental 
compositions similar to the Mars soil, only iron smectites, and particularly 
montmorillonite, reasonably simulated the chemical reactivity in the Viking 
Biology LR experiment. These simulations showed that if any CO2 was to be 
released in the LR Viking experiment, calcite (CaCOg) could not be present in 
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the dusty surface materials at concentrations higher than about 0.5%, and the 
pH of the soil could not be above ca. 5. 

In separate simulation experiments, Hubbard (16), using the modified 
clays of Banin, detected photochemical fixation of atmospheric and its 

incorporation into organics in the soil. This photochemical fixation was 
similar to the observations on Mars in the Viking Pyrolytic Release (PR) 
experiment. Other candidates for Mars soil components which were suggested 
using Viking Biology results as the guideline were carbon suboxide (17) and 
various mafic silicate minerals (18). 

Spectroscopic evidence: Earth-based visible and near infrared 

reflectance data were accumulated through telescopic measurements of McCord's 
group in Hawaii (19,20). Since the termination of the Viking extended mission 
in 1978, this spectral information became the major if not the sole source of 
new experimental data on Mars soils mineralogy. The Mars reflectance spectrum 
was found to show opacity and band saturation (no reflectance) in the near UV- 
visible range (0.3-0. 5 ym ) with increasing transparency (increasing 
reflectance) in the long-wavelength visible and near infrared range (0.5-0. 9 
pm). However, it was lacking any clear and pronounced absorption features 
throughout this range. Spectra of crystalline, pure iron oxides such as 
hematite and goethite were found by Singer (21) and by Sherman et al. (22) to 
exhibit significant deviations from the Mars reflectance spectrum measured 
from Earth. However, several iron-containing minerals have been found to bear 
significant similarity to the Mars reflectance spectrum. These include 
amorphous iron-silica gels (23,24), palagonite (21), iron-containing aluminum 
oxy-hydroxide (25), and iron saturated montmorillonite clay (14). A feature 
common to all of the candidates is that they contain ferric iron in a matrix 
of oxy-hydroxide which is not well crystallized, or has only short-range 
ordering. Unfortunately, on the basis of the spectral comparison alone, there 
is yet no obvious way to choose among the various candidates, and additional 
properties and characteristics of the Martian soil and its Earth analogs have 
to be studied in order to further constrain the possible Martian mineralogical 
composition. 
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IMPLICATIONS OF A POLAR WARMING FOR DUST TRANSPORT INTO 
THE NORTH POLAR REGION OP MARS; J.R. Barnes and J.L. Hollingsworth, 
Department of Atmospheric Sciences, Oregon State University, Corvallis, 

OR 97331 

The layered terrains in the polar regions of Mars have been widely in- 
terpreted as records of quasi -per iodic variations in the Mars climate system 
driven by orbital element variations (see, e.g.. Pollack and Toon, 1982). 

It has been suggested that the polar leiminae in the north are being formed 
at present and that perhaps the bulk of the deposited dust and water ice is 
transported to the pole during global dust storms (Cutts, 1973; Pollack et 
al., 1979), which currently occur during the growth phase of the seasonal 
CC^ frost cap. Pollack et al. hypothesized that CO 2 and water ice condensa- 
tion onto dust grains would consitute a very effective mechanism for remov- 
ing the dust to the surface in the polar regions. Paige (1985) carried out 
analyses of Viking IRTM observations and concluded that deposition of dust 
did occur (probably in conjunction with (X >2 condensation in the atmosphere) 
in the region of the north residual cap during the winter solstice (the sec- 
ond) dust storm of 1977. Die magnitude of this deposition is uncertain, 
however, Jakosky (1983) and Jakosky emd Martin (1984) have argued that the 
transport of dust (and water) onto the north residual cap during global dust 
storms is insignificant. 

The deposition of dust onto the seasonal polar cap is also of consider- 
able importance, because of the effects on the polar cap albedo and thus the 
radiation budget of the subliming cap. Additionally, the dust loading of 
the polar atmosphere is significant because of its effects on the radiation 
budget of the cap: the enhanced emissivity of the atmosphere implies en- 
hanced downward IR fluxes. 

The intense warming of the polar atmosphere which was observed by the 
Viking IRIM during the winter solstice dust storm of 1977 (Martin and Kief- 
fer, 1979; Jakosky and Martin, 1984) must have been the result of substan- 
tial dynamical transports of heat into the north polar region. The warming 
at the pole (-50 K in magnitude in the vicinity of the 25 km level) must 
have been produced by congressional heating associated with downward - in a 
Lagrangian sense - motion of air parcels. An accompanying poleward flow is 
then also implied. Such a Lagrangian motion pattern implies poleward and 
downward transport of tracers such as dust emd water. 

Numerical simulations conducted with a simplified beta-plane wave-mean 
flow model strongly suggest that the polar warming could have been induced 
by a forced, vertically propagating, planetary-scale (zonal wavenumber 1) 
wave (Barnes and Hollingsworth, 1985). The dynamics of such a warming are 
roost clearly viewed in the so-called transformed Eulerian-mean system, in- 
stead of the traditional Euler ian one. In the transformed system, the warm- 
ing is produced by downward vertical motions - rather than eddy transport of 
heat. [The transformed Eulerian-mean vertical motion field is generally a 
good approximation to the true Lagrangian-mean vertical motion field - un- 
like the Eulerian-mean vericzd. motion field.] The mean meridional circula- 
tion in the transformed Eulerian-mean systoa, the so-called residual mean 
meridional circulation, is characterized by pole%«ard and downward flow into 
the polar interior during the simulated warmings. This poleward and down- 
ward flow takes place throughout a very deep region as the warming occurs, 
with the much weaker return flow confined to low levels. The residual mean 
meridional circulation should constitute a good first approximation 
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appropriate to the transport of tracers (Holton, 1981), thus implying 
transport into the polar interior of a tracer initially located largely 
equatorward of the cap boundary (e.g., dust or water). 

Pull tracer (dust) transport simulations corresponding to the dynam - 
ical warming simulations are planned. These will allow a quantitative 
examination of the transport of dust into the polar region in association 
with a polar warming event. 
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DUST DEPOSITION AND EROSION ON MARS: CYCLIC 

DEVELOPMENT OF REGIONAL DEPOSITS 

Philip R. Christensen, Dept, of Geology, Arizona State University, Tempe, AZ, 85287. 

The cycle of dust erosion, transport, and deposition plays an important role in the evolution 
of the martian surface, both today and in the past. Thermal, radar, and visual remote sensing 
observations provide important constraints on the surface properties, and have been used to 
determine the location and physical properties of regional dust deposits. These deposits provide 
direct clues to the rate and history of dust deposition, and suggest that the martian surface is being 
actively reworked. The observed dust cycle may be directly related to the cycle of variations in 
Mars' orbit, with dust continually transported between hemispheres on time scales of 10^ to 10^ 
years. 


Major dust deposits are currently located in three northern equatorial regions: Tharsis 
(-20°S to 59°N, 600 to 190°W), Arabia (-S^S to 30 pN, 300° to 360oW), and Elysium (lO® to 
30°N, 21QO to 225°W). They are covered by fine (~2-40 p.m), bright (albedo > 0.27) particles, 
with fewer exposed rocks and coarse deposits than found elsewhere (1,2). These regions also 
have less bonded material exposed at the surface than found elsewhere (3). Dust deposits may 
have initially formed due to differences in wind velocity. Once initiated, the burial of sand and 
rocks would make removal of fines increasingly difficult, enhancing the rate of dust accumulation. 
Dust is currently deposited uniformly throughout the equatorial region at a rate of ~40 |Xm/global 
storm. Over geologic time the rate of accumulation may vary from 0 to 250 |i.m/year due to 
changes in atmospheric conditions produced by orbital variations (4). Dust deposited during 
glob^ storms is subsequently removed only from dark regions, resulting a net accumulation in the 
low-inertia, bright regions. The evidence for the subsequent removal from dark regions comes 
from: 1) the observed higher dust content over dark regions than bright regions during clear 
periods; 2) the post-storm darkening of dark regions (5, this issue); 3) the removal of storm- 
deposited dust at the Viking lander 1 site during the year following the major 1977 dust storms (6); 
and 4) the historical persistence of classic dark features. Non-removal of dust from low-inertia, 
bright regions results in a net accumulation of dust in these areas. 

The thickness of the current dust deposits can be estimated from thermal, radar, and visual 
observations.. The low thermal inertia of the deposits places a lower limit of ~0.1 m on their 
thickness, while the sparse but ubiquitous presence of exposed rocks and the degree of visible 
mantling indicates that the thickness is less than 5 meters (7). Dual-polarization radar observations 
of Tharsis reveal a very rough texture. These measurements can best be reconciled with the other 
observations by assuming that a relatively thin dust layer buries most of the surface rocks but that 
the radar samples through this layer to the rough surface below (1). The maximum thickness of 
this layer can be estimated from the electrical properties of rock powders. For dry powders, the 
radar energy should penetrate 5 to 50 wavelengths (0.6 to 6 m for the 12 cm radar used) before 
being attenuated by 1/e (8). The presence of adsorbed water will not affect these results (9). 
Thus, 12 cm radar could penetrate ~ 1 m, be reflected, and exit the surface with only a 1/e 
reduction in initial energy. This estimate puts a constraint of 1-2 m on the dust mantle thickness, 
which agrees well with the value obtained using the exposed rock abundance derived from thermal 
measurements. 

Based on their thickness and rate of accumulation, the age of these deposits is 10^-10^ 
years, suggesting a cyclic process of deposition and removal (1). One possible cause may be 
cyclic variations in the magnitude and location of maximum wind velocities related to variations in 
Mars’ orbit. At present, perihelion and maximum wind velocities occur in the south whereas 
regional dust deposits occur in the north, suggesting net transport from south to north. Orbital 
parameters oscillate with periods ranging from 5x10^ to 10® years. The agreement between these 
periods and the dust deposit age suggests that there is a possible link. At different stages in orbit 
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evolution, maximum wind velocities will occur in the north, with subsequent erosion and 
redistribution of the accumulated fines. 

The model proposed here implies that material must be periodically removed from 
regional dust deposits in order to prevent long-term buildup of fine material in a given location. 
Thus, the proposed young age of these dust deposits requires mechanisms for eroding extensive 
deposits of fine material, llie burial of sand and rocks makes it increasingly difficult to set 
particles in motion, with 20 ^tm particles requiring wind velocities a factor of 2 higher than the 
most readily moved particles (10). There are, however, several mechanisms for eroding these 
deposits, including erosion from the edge inward, and increased surface shear stress produced by 
increased winds (1). 

A more plausable mechanism involves the formation of coarse particles as bonded 
aggregates of dust. Bonding of material has been observed at both Viking landing sites, and 
globally pervasive crusts have been detected from remote sensing observations (2). These crusts 
may form during transition periods between obliquity extremes as volatiles carrying adsorbed ions 
are cycled into and out to the surface (2). This mechanism would link the erosion of the deposits 
to the same process that leads to their formation. Thus, crusts could form after the deposition of a 
layer of fine material has been deposited after each obliquity cycle. Such crusted aggregates could 
provide a source of coarse particles that could be more readily moved by the wind, thereby 
providing a mechanism for eroding the underlying deposit of dust 

Perhaps the most important process may be surface erosion due to insolation-driven 
convective vortices (dust devils) of various scales. Experimental work suggests that vortices are 
very effective in raising particles of all sizes which can then be easily transported by much lower 
winds (11). There is direct evidence for the occurrence of dust devils on Mars, both through the 
passage of 5-950 m diameter vortices at the Viking lander sites (12), and through direct 
observation of dust devils from orbiter images (13). Of the 118 vortices observed at die lander 
sites, 4 had wind velocities greater than 30 m/sec, which may have been sufficient to raise dust 
(12). These vortices product a factor of 2-3 enhancement in die ambient wind velocity. Because 
dust devils form due to convection in an atmosphere with a superadiabatic lapse rate, they are more 
frequent during periods of maximum surface heating. This erosion mechanism would be most 
effective over low-inertia surfaces during summer, and may provide a mechanism for eroding dust 
from the hemisphere that has the maximum solar insolation during summer. Thus, dust deposited 
in one hemisphere when the insolation maximum was in the opposite hemisphere would 
subsequendy be eroded when obliquity variations caused the insolation maximum to reverse. 

The presence and history of regional dust deposits provide evidence for cyclic processes of 
deposition outside the polar regions and support models of cyclic variations in martian climate over 
geologic time. Dust is continually eroded and redeposited, with the location of major deposits 
shifting on time scales of 10^ to 10^ years. In this model much of the uppermost martian surface is 
young and is being constandy reworked. 
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SURFACE ALBEDO VARIATIONS ON MARS: IMPLICATIONS FOR 
YEARLY DUST DEPOSITION AND REMOVAL 

Philip R. Christensen, Dept of Geology, Arizona State University, Tempe, AZ, 85287 

Dust deposition and removal is an important process in the development and modification 
of the martian surface. Mars has been known to have variable surface markings from the earliest 
telescopic views of the planet These changes have since been seen to be related to aeolian activity, 
primarily through tlie reworking of bright dust deposited following major global dust storms (1,2). 
Viking Infrared Thermal Mapper (IRTM) observations of albedo have also revealed significant 
changes in surface brightness through time, again primarily associated with major global dust 
storms (3,4). All •' f se observations indicate that there is a significant amount of dust that is 
deposited during the decay of global storms which is subsequently reworked and redistributed. 
The purpose of this study is to determine the degree, spatial distribution, and timing of the 
deposition and removal of dust-storm fallout, and to relate the current patterns of dust deposition 
and removal to the long-term evolution of the martian surface. 

A model has been proposed (5, 6 this issue) for the development of regional dust deposits 
that form through the preferential accumulation of dust-storm fallout into specific northern 
hemisphere regions. In this model, dust is deposited uniformly during the decay phase of each 
major storm, but is subsequently removed only from regions that are seen today as classic dark 
areas. Thus, dark regions remain unmantled by dust, whereas bright regions have developed a 1-2 
m thick mantle of fine, bright dust (5). This model can account for the Wgh thermal inertia (coarse) 
material observed in dark regions, together with their relatively high rock abundance (7), and low 
albedo. Conversely, bright regions have fine particles (5-40 pm) and fewer exposed rocks, 
presumably due to mantling of the coarse material by dust 

In order to directly observe the seasonal changes in surface brightness associated with dust 
deposition and removal, the albedo of specific regions in both hemispheres has been determined 
through time. The IRTM data were collected into 1° latitude by 4° longitude bins, at 3 hour 
intervals for each 10° of Lj. Using these data, the albedo changes for a given area have been 
investigated from the beginning of the Viking mission (Lg 84°), through the first (Lg 190-240°) 
and second (L, 270-340°) global dust storms that occurred in 1977. Global data are available 
through Lg 120° of the second year, allowing a year to year comparison of surface albedo. 

The albedo variations as a function of season are shown in Figure 1 for representative 
bright and dark regions. All of the areas studied show a marked increase in brightness associated 
wife fee two global storms, due primarily to fee presence of dust in fee atmosphere. The increase 
in brightness, even for bright regions, indicates feat fee albedo and scattering phase function of 
suspended dust varies from dust on the surface. The maximum brightness at the peak of fee 
second storm was nearly equal for most bright and dark regions, indicating feat fee atmospheric 
dust was optically thick. For some dark regions, however, such as Solis Planum, fee albedo 
remained relatively low even at fee height of fee storm activity, suggesting feat fee atmospheric 
dust was not globally uniform nor well mixed. Many areas show a non-uniform decrease in 
brightness during fee decay phase, again suggesting spatial variations in dust load and non- 
uniform mixing, possibly due to episodic injection of dust into the atmosphere locally (8). 

The albedo of most regions had returned to the pre-storm value by Lg 355°, indicating feat 
fee atmosphere had cleared to pre-storm levels by feat time. This conclusion is supported by 
Viking lander observations, which show feat fee opacity over fee two lander sites had decreased to 
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pre-storm levels by Lj 360® (8). Therefore, surfaces that remained brighter after 360° than they 
were prior to the two storms are thought to be covered by a thin layer of bright dust fallout. 

The distribution of surfaces that remained bright following the storms, and those where the 
surface quickly returned to its pre-storm albedo follow a consistent pattern. The albedo of bright 
regions, such as Arabia and Tharsis, rapidly returned to pre-storm values, and was close to the 
albedo of the previous year (Fig. la). Many dark regions also darkened to nearly their pre-storm 
levels by Lg 360° (Fig. lb). This pattern holds particularly well for southern hemisphere dark 
regions. This behavior is consistent with the model of deposition described above; in dark regions 
the dust is rapidly removed with little net accumulation, whereas in bright regions a dust mantle 
already exists so that the deposition of additional bright dust does not affect the surface albedo. 

There are several dark regions that differ from the general trends described above and 
provide insight into the level of dust activity that occurs throughout the year. Syrtis Major and 
Acidalia Planitia are among the few regions that remained significantly brighter at Lg 360° than they 
were before the global storms began. These areas did, however, continue to darken with time, 
returning to nearly their pre-storm albedo by Lg 120° (Fig. Ic). It is interesting to note that the 
albedo of these and some other regions was still slightly higher at this time than it was the previous 
year, suggesting that some dust still remained on the surface. This finding is consistent with 
observations at the Viking lander 1 site where dust was deposited following the global storms and 
was not removed until over a year later (9). These observations support the hypothesis that Syrtis 
Major and Acidalia Planitia act as local dust sources during inter-storm periods, producing 
enhanced dust loading in the northern hemisphere (10). 

In summary, observations of seasonal changes in surface albedo reveal regional differences 
in the deposition and subsequent erosion of dust-storm fallout. Southern hemisphere dark areas 
quickly return to close to their pre-storm albedos, suggesting rapid removal of any dust that was 
deposited. Northern hemisphere dark regions are brighter post-storm, but gradually darken to pre- 
storm levels over a Mars year. In doing so they act as locd sources of dust during otherwise clear 
periods. Dust does not appear to be removed from bright regions, resulting in the 1-2 m thick 
deposits observed today. 
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1.2.1 FOURIER COLLOCATION 

For problems with periodic boundary conditions, the Fourier expansion of 
I a function u(x) is given by the infinite series 


u(x) = I a^e 

k=-oo 


( 8 ) 


The collocation projection is defined by the discrete Fourier transform pair 


N/2-1 . 

= I 

^ k=-N/2 ^ 


ikx 


(9) 


where the coefficients are defined by 


1 N-1 -ikx . „ „ 

^k ^(Xj)e k 2* 2"^^*****2'”^' 


( 10 ) 


The collocation points, x ^ , are uniform on the interval [0,2ir] 


Xj = 2irj/N j = 0,1,2,...N-1. 


( 11 ) 


The transforms (9) and (10) are almost always computed by the use of a fast 
Fourier transform if N is a highly composite integer such as N = 2^39. 

Derivatives of the function u at the collocation points are approxi- 
mated by the derivatives of the interpolating polynomial. Thus, the Ith 
derivative of u is approximated 


d*'P„u N/2-1 

— l 

dx* k— N/2 


xl,'' ikx 
(Ik) u,,e 


( 12 ) 
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From the form of equation (12), it is clear that the evaluation of the deriva- 
tive at the collocation points can also be computed efficiently with a fast 
Fourier transform. See Hussaini, Streett, and Zang [4] for more information 
on the implementation of the Fourier collocation method. 


1.2.2 CHEBYSHEV COLLOCATION 

The collocation points for using a Chebyshev method to approximate a non- 
periodic function are usually defined by 


Xj = -cos(irj/N) j = 0,1,..., N. 


(13) 


These points are the extrema of the order Chebyshev polynomial, T^(x), and 
are obtained from the Gauss-Lobatto Integration formula (see Davis and 
Rabbinowitz [5]). 

The collocation projection operator is defined as the interpolation 

V "„Vn« 

n~0 


where the coefficients are defined by 
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where 


(2 j = 1,N 
( 1 otherwise 


(15) 


Again, derivatives of u at the collocation points are approximated by 
the derivative of the interpolating polynomial evaluated at the collocation 
points. The first derivative, for example, is defined by 


dx n 
n=0 


I„(x) 
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( 1 ) 

N+1 


= 0 
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n ti n+2 
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( 16 ) 


The transform pair given by equations (14) or (16) and (15) can be 
efficiently computed with a fast cosine transform. Equivalently, the inter- 
polating polynomial and its derivatives can be computed using matrix multipli- 
cation. The matrices for the Chebyshev collocation method are conveniently 
collected in the review by Gottlieb, Hussaini, and Orszag [6]. For N < 32, 
this approach is competitive with using a fast cosine transform, at least on 
serial computers. 


1.3 APPROXIMATION THEORY (COLLOCATION) 

1.3.1 FOURIER COLLOCATION 

The problem of how well Pjju approximates u for Fourier approximations 
has been discussed by Kreiss and Oliger [7], Pasciak [8], and by Canuto and 
Quarteroni [9]. See also Mercier [10]. It is most convenient to express the 
interpolation results in terms of a Sobolev space, H™(0,2ir). This is a 
Hilbert space with the norm 
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lul - I |u|^ (17) 

^ j=0 J 

defined in terms of the seminorms 
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(18) 


The use of the discrete Fourier transform pair (9), (10) represents the 
projection of the Sobolev space onto the space Sjj(0,2n), the space of 
Fourier polynomials of degree N. 

The primary interpolation result is given by Theorem 1: 


Theorem 1: For any 0 ^ p ^ q with q > 1/2 there exists a 

constant C independent of u and N such that 


Hu - P„ull < CnP"‘1|u| . 

N p — q 


(19) 


Proof ; See Pasciak [8]. 

Equation (19) states that the rate of convergence depends (through the 
order of Sobolev norms) only on the smoothness of the function being approxi- 
mated. This type of error decay is known as spectral accuracy. In practice, 
one sees errors which decay exponentially and hence spectral accuracy is often 
called exponential accuracy. Several applications described in Section 2 
exhibit exponential accuracy. The term infinite order accuracy is also used 
often to refer to the case as q 
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Exponential accuracy has been shown explicitly by Tadmor [11] for func- 
tions u which are also analytic in the complex plane* 

Theorem 2: Let u(x) 2Tr-periodic and analytic in a strip of 

width 2sq. Then for any 0 < s < sq 

llu - P„ull < CM(s)/sinh(s)N^e (20) 

N p — 

where C depends on p and 


M(s) = Max |u(z) I . 

I Imz I £ s 

Proof : See Tadmor [11]* 

If the solution is not very smooth, then the approximation may not be 
very good. In fact, if the function is discontinuous, the interpolant shows 
global oscillations (Gibbs phenomenon) and the approximation error decay is 
globally only first order. Smoothness is not usually a problem with the solu- 
tions of many elliptic or parabolic equations, but discontinuities are 
characteristic of the solutions of hyperbolic equations. 

It is still possible to obtain spectrally accurate approximations to non- 
smooth functions, at least away from any discontinuities, but some type of 
filtering is required. Two papers which address this issue are Majda, 
MacDonough, and Osher [12] and Gottlieb and Tadmor [13]. The first approach 
used to smooth discontinuous solutions was that of Majda, MacDonough and Osher 
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[12] whose results show that spectral accuracy can be retained if Fourier 
space filtering is applied. Since the main results refer directly to the 
solutions of hyperbolic partial differential equations, they will be discussed 
in the next subsection. 

Gottlieb and Tadmor [13] have taken the approach of smoothing in real 
space to allow the accuracy to depend on the local smoothness of the func- 
tion. The smoothing procedure consists of convoluting the collocation approx- 
imation with a regularization kernel which is localized in space. If we 
call P u the smoothed approximation to the originally oscillatory inter- 

polant Pfju, the convolution takes on the form 


Pu(x) 



N-1 


I 

j=0 


PjjU(yj)i|»®’P(x - yj) 


( 21 ) 


where 


<l»®’P(y) = 


1 ^y-. 


sin((p + -|>|-) 
sin(y/6) 


( 22 ) 


is the Dirichlet kernel localized in space by the cutoff function 


^a5^/(5^-l) 

0 


| 5 | < 1 

otherwise 


(23) 


The function p ensures that the kernel does not Interact with any regions 
of discontinuity* For example, for a single discontinuity at x = it, they 
choose 0 = ir|x - tt|* With this smoothing, they show that the error depends 
only on the smoothness of the cutoff function p(?): 
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Theorem: Let p (5 ) be a cutoff function satisfying 

p(0) = 1 and having support In [”TT,Tt]. Then for any x [0,2ir] 

the smoothed function Pu satisfies the estimate 


|Pu(x) - u(x)| < C llpll Max|D’*^u(y) |(1 + 0 s > 1. (24) 

” ® |y-x| <6Tr 

0 < k < 2s 


Proof; See Gottlieb and Tadmor [13]. 


1.3.2 CHEBYSHKV COLLOCATION 

To study the approximation properties of the Chebyshev projection (14), 
it is practical to work in a weighted Sobolev space with weight w(x) = (1 - 
x2)”l/2. Defining the weighted L^ norm by 


II ull 


2 

0,w 


(u,u) 


w 


= j u^wdx 
-1 


(25) 


and the Sobolev norm by 


II uH 


2 

q,w 


3 .A. 2 


I 


i=l dx 


r 0,w 


(26) 


the spectral approximation result is given by 


Theorem 4: Let q >1/2 and 0 ^ p ^ q. 

constant C such that for all u in H^(— 1,1) 


Then there exists a 
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lu - P„ull < CN 
N p,w — 


2p-q 


lull 


q,w 


(27) 


Proof ; See Canuto and Quarteroni [9]. 

So, like the Fourier approximation, the Chebyshev interpolation gives 
spectral accuracy; that is, the accuracy depends only on the smoothness of the 
function to be interpolated. Exponential convergence has also been proved by 
Tadraor [11]. This time, the function u must be analytic in an ellipse with 
foci at -1 and 1: 

Theorem 5: Assume u(x) is analytic in [-1,1] and has a regularity 

ellipse whose sum of its semi-axes equals r^ - exp(riQ) >1. Then for any 
n,0 < n < tIq we have 

1/2 

llu(x) - P ull , < Ne~^ (28) 

- 1 

where the norm is defined by 


Hun . = ^ (n-p)^lu 1^. 

p=0 P 


Proof ; See Tadmor [11]. 


(29) 


If the function which is being approximated is discontinuous, it is still 
theoretically possible to recover a spectrally accurate solution [13] by fil- 
tering in physical space. The procedure is the same as the smoothing proce- 
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dure for the Fourier case, but the Dlrlchlet kernel is replaced by 


Kp(y) 


. o - 

PIT 


i;(y) 

“T~ * 


( 30 ) 


1.4 THEORY OF SPECTRAL COLLOCATIOH METHODS FOR PDE'S 

Proofs of the convergence of spectral approximations to partial differen- 
tial equations are usually accomplished using energy methods which mimic 
proofs of the well-posedness of the original equations. Consequently, it is 
most convenient to discuss stability and convergence with respect to the three 
major types of partial differential equations separately. 


1.4.1 ELLIPTIC EQUATIONS 

Theoretical analysis of the convergence of Fourier collocation methods is 
simplified because of periodic boundary conditions. The elliptic problem is 
to find the function u(x) which satisfies 

Lu = f x€[0,2ir] (31) 

u(0)=u(2tt ) 


where L has the property 


(Lu,u) > a Hull 


2 

1 


a > 0. 


(32) 


The Fourier collocation approximation is obtained as described in section 
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1.2.1 and satisfies the same inequality, i.e., ¥ueSjg 


(L u,u) > allull 
c — 1 

Then we have the 


Theorea 6: For t > 1 

f€H^"^(0,2it) and u€H^(0,2ti) 

P • P ’ 


there exists a constant C such that If 
then the following estimate Is optimal; 


lu - < C(1 

N 1 — 


+ n2)0-t)/2||ui| 


(33) 


Proof ; See Mercier [10]. 


Chebyshev methods with both Dirichlet and Neumann boundary conditions 
have been analyzed for the elliptic differential equation of the form 


Lu = -(au ) + (bu) • 

XX X 


(34) 


The Chebyshev spectral collocation approximation is formally written as 


Vc = - ^ 


aPN(bu^) 

+ — 35 ; — + 


(35) 


For Dirichlet problems, the equation is collocated at the interior points and 
boundary conditions of the form 


u(-l) = u, and u(+l) = u 
1 r 


( 36 ) 
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are specified directly at the boundary points. Stability and convergence were 
proved by Canuto and Quarteroni [1^1 using a variational approach. They show 


7* bet Ug be the solution to LgU^. = where L^. is defined 
by equation (35) with homogeneous boundary conditions, u^, = uj = 0 then with 
suitable conditions on a, b, o the following estimate holds 


llu - u II, < C,N^“’^llull + C-N"®«fll 
c 1 ,w — 1 r,w 2 


s,w 


(37) 


Proof ; See Canuto and Quarteroni [14], Theorem 2.4. 

Convergence proofs for Neumann or mixed-type boundary conditions are 
available for boundary conditions applied in one of two different ways. A 
discussion of these approaches can be found in Canuto [15], [16]. The first 
approach is explicit. At interior points, the equation is collocated normally 
as in equation (35). At the boundary points, however, the collocation approx- 
imation to the derivative is written in matrix form and the boundary condi- 
tions are used to determine the value at the boundary point. Thus, the 
approximation to the boundary condition 

B,u = 3u + au 
1 X 

(38) 

B u = 6u + yu 
r X 


is found by solving the system 
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N-1 




■'■ '*' ‘^NO^O 


( 39 ) 


where Uj = u^.(xj) and is the matrix for the derivative at the 

collocation points (see Gottlieb, Hussaini, and Orszag [6]). 

The convergence is very rapid for smooth solutions: 


Theorem 8: Let a > 1/2 and let u and be solutions to Lu = 

f and L^u^ = f^ where L and are defined as above > Then with 

explicitly applied Neumann boundary conditions the following convergence 
estimate holds 




(40) 


where n = (1 ~ x'^)w(x) and C is independent of N« 


Proof : See Canute and Quarteroni [14], Theorem 3.2* 


Canute [16] also describes how to impose Neumann boundary conditions 
implicitly for elliptic problems. In this way, the boundary conditions are 
not exactly satisfied because what is actually solved is the modification of 
the interior approximation. For the spectral case of a pure Neumann problem, 
the first derivatives are computed normally as in equation (16). At the 
boundary points, the derivatives are replaced by the Neumann conditions. Then 
the second spectral derivatives are computed by using (16) again on the modi- 
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fied set of derivatives. This has the advantage that all of the points are 
treated the same, but the boundary conditions are not exactly satisfied. The 
boundary error does decay spectrally, however. 

Theorem 9: Let u^ be the solution to L^u^. = f with implicit Neumann 

boundary conditions. If u€H™(-l,l) with m > 5/2 then 

1^ (±1)| < CN^“"llull (41) 

'3x ' — m,w 

where C > 0 is Independent of N. 

Proof: See Canute [16]. 

The convergence in the interior is also spectral, and the estimate bounds 
both the solution and the collocation derivative. 


Iheoren 10: Under the assumptions of Theorem 9, 


N-1 

f 

j 


V0,w-^ 1 


" CN^"™Hull (42) 

m,w 


where the Wj are the Gauss-Lobatto weights at the points x j 


Proof: See Canute [16] 
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1.4.2 PARABOLIC EQUATIONS 

The convergence and stability theory for linear parabolic equations, like 
the theory for elliptic equations, is fairly well developed. In particular, 
the theory has centered on studies of semi-discrete equations in which the 
spatial variation is discretized, but the time variation is left continuous. 
The emphasis, however, has been on application to boundary value problems — 
that is, on the convergence of Chebyshev collocation methods. In this section 
we survey theoretical results for initial-value problems of the form 


u^ = (Au ) + Bu + Cu + f 

t XX X 

(43) 


u(x,0) = 


where A, B, and C are n x n matries. The general collocation approxima- 
tion to the first, third, and fourth terms of the right hand side of equation 
(43) is written in a manner identical to that of the elliptic equations in 
equation (35). 

Stability of the Fourier approximation of the heat equation is easy to 
prove and is discussed in Gottlieb, Hussainl, and Orszag [6]. The more com- 
plicated case is equation (43) above. Kreiss and Ollger [17] have proved 
stability with two different treatments of the first order term. The first 
treats it in "skew symmetric form", that is, by writing 

Bu + Cu = i (Bu + (Bu) ) + (C - i B )u. (44) 
X 2 X X 2 X 


The second, discussed more fully in the next section, involves filtering the 
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flrst derivative to ensure stability. A convergence estimate using the skew 
symmetric form for the scalar equation with f = 0 is 

Theoren 11: Let t 1, T > 0, and assume that UQeHp^^(0,2ir) . Then 

there exists a constant C such that the following estimate holds: 

llu(t) - u^IIq < C(1 + [IIUqIiJ + (HUqII^+ 

for all t£[0,T]. 

Proof : See Mercier [10] Theorem 11.2. Here Hp is defined in terms of 

distribution derivatives of periodic functions. 

The convergence of Chebyshev approximations to parabolic equations on 
bounded domains has received a lot of attention recently. The spatial approx- 
imation for a Dirichlet problem will be exactly like that for the elliptic 
problem. Stability for the heat equation with non-constant coefficients was 
originally shown by Gottlieb [18]. Convergence estimates were worked out by 
Canuto and Quarteroni [19]. For the scalar heat equation 

u^ = a(x)u X e (-1,1) (46) 

t XX 

with homogeneous boundary conditions u(-l,t) = u(l,t) = 0 they show 

Theorem 12: Let a > 1/2 and S > a +2 and T > 0. If 
1 s 

u£L (0,T; H (-1,1)) then there is a constant C, independent of N such 

W I...- — 



-2 



Neumann boundary conditions for parabolic problems can also be applied 
either explicitly or implicitly. For the implicit treatment, convergence is 
similar to that of the corresponding elliptic equation: 


Theorem 13: Suppose the solution to the differential equation (46) with 

Neumann boundary conditions is regular to the extent that u€L^(0,T;H^(-l , 1 ) ) 
and the time derivative satisfies u^€L^(0,T;H^^(-l,l)) for T > 0 and m 
>5/2. If_ Uq€H®(- 1,1) then 
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1.4.3 HYPERBOLIC EQUATIONS 

The study of the convergence of spectral approximations to hyperbolic 
equations is complicated by the fact that the straight-forward discretization 
of an equation of the form 


written as 


u^ + a(x)u + bu = 0 

t X 


d u du 

FT + bu^ = 0 


(50) 


(51) 


is often unstable. In this section, we will discuss the available theory of 
formally stable approximations. 

Fourier methods are stable if a(x) is of fixed sign. If a(x) in 
equation (50) is strictly positive and b = 0, the energy estimate for the 
approximation , (51) 


d . 0 


(52) 


shows that the approximation is stable. If a(x) is zero at some point, how- 
ever, then this estimate is not valid and no general technique is available to 
show stability. 

Two basic approaches have been used to devise schemes which can be shown 
to be stable. The first, indicated in the last section, is to write the spa- 
tial derivative in skew-symmetric form. That is. Instead of computing (51), 
one computes 


Tt~ 


1/2 


ad u 


c ^ 5(PNaUc> 


- 1/2 P,{^} - bu^ = 0 


(53) 
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Kreiss and Oliger [20], [17] showed that this discretization is stable. 
Mercier [10] examined the stability and convergence of the Fourier approxima- 
tion to the skew-symmetric equation 

u^ + v(x)u^ + (v(x)u)^ = 0 (54) 

and showed that the error decay is spectral. 

Theorem 14: For t > 1 and T > 0, if the initial condition 

satisfies u(x,0)€Hp^ (0,2ir ) then there is a constant C independent of N 
such that 

llu(t) - u (t)ll < C(1 + " '^^/^llu.ll . (55) 

C — U T 


Proof ; See Mercier [10], Theorem 9.1. 

Though approximations written in skew symmetric form are stable, there 
are objections to their use. The first objection is that they are less 
efficient since they have twice as many derivatives to evaluate. More 
important, conservation is lost when this is applied to conservation law equa- 
tions (such as the equations of gas-dynamics) for the computation of weak 
solutions. Tadmor [21] has examined the skew-self adjoint form of systems of 
non-linear conservation laws. They can be explicitly shown to be well-posed, 
but the conservation property is lost. 

The alternative to rewriting the equation in skew-symmetric form is to 
use the approximation of equation (51) and filter the solutions. Finite 
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difference solutions are often filtered by adding an explicit low order 
artificial viscosity. The goal of filtering Fourier spectral solutions is to 
do so without destroying the accuracy of the method. 

Two approaches for filtering Fourier approximations to guarantee 
stability have been suggested. The first was proposed by Majda, McDonough and 
Osher [12]. In their method, the spectral derivative defined in equation (12) 
is modified by filtering the computed solution. For linear problems, this can 
be done efficiently by modifying the Fourier coefficients of the solution and 

using those new coefficients when the derivative is computed. Let 

00 

p(x)€c (-ir,x) be a "filter function". Its values are zero near 9 = ± ir 
and identically one in a neighborhood of 6=0. The Fourier coefficient 

^ /s 

u^ is replaced by p(2irk/N)u^ and this is used in equation (12) to 
compute the derivative. 

For smooth initial conditions, smoothing gives a stable approximation and 
spectral accuracy 

Theoroi 15: For the error satisfies the inequality 

llu(x,t) - u.(x,t)y. < C h^ for all s,X (56) 

where C depends on both s and X. 

Proof : See Majda, McDonough, and Osher [12], Corollary 1. 


If the solution is discontinuous. It is still possible to obtain spectral 
accuracy in the sense of equation (56) if the initial condition is properly 
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smoothed. It is not enough to smooth the discrete Fourier coefficients of the 
initial condition with a filter whose support is enclosed within the support 
of p. Rather, it is necessary to use smoothed versions of the exact 
Fourier coefficients. 

A different approach to filtering was proposed by Kreiss and Oliger 
[17], Instead of filtering the solution with a predefined filter, they showed 
that linear stability could be obtained by smoothing the space derivative with 
a weak filter which depends on the smoothness of the function. They arbi- 
trarily split the frequency range of the solution into a high frequency range, 
|k| > Nj, and a low frequency range, |k| < Nj, The coefficients of the low 
frequency range are not modified at all. The coefficients of the high fre- 
quency range are modified only if they do not decay rapidly enough. Call 
v(x) = defined by equation (12) and define vj to be the derivative 
summing only the low frequency components |k| ^ N^. The modified 
coefficients of the derivative are defined to be w = Hv where 


”k M ^ 


for |k| £ Nj^ 

if |k| > and |v^| < D Hv^ll/|k|^ (57) 


D «VjllVj^/(|vj^| |k|^) 


otherwise. 


Kreiss and Oliger prove the following stability theorem: 


Theorem 16: Suppose the coefficient a(x) in equation (50) is smooth so 
that its Fourier coefficients decay at a rate |k| The approximation 


H3u- 

“t * -5T-> = 0 


( 58 ) 
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where the filter H is defined by (57) Is stable and converges If j = 3 > 2. 

Proof ; See Kreiss and Oliger [17], Theorem 4.2. 

For linear problems, however, it is not clear that filtering is always 
needed. The fact that the energy method gives only a sufficient condition for 
stability means that equation (52) does not prove instability if a is not of 
one sign. For example, Gottlieb, Orszag, and Turkel [22] show stability in 
the usual sense of convergence as N -»• «> of the scalar equation where a(x) 
= Asin(x) + Bcos(x) + C for arbitrary A, B, C. The numerical solutions do, 
however, grow in time - just as the exact solutions do. 

For non-linear problems, experience shows that filtering of the Fourier 
approximation is needed, particularly if there are discontinuities in the 
solutions. Hussainl, Kopriva, Salas, and Zang [51] discuss the application of 
these filtering methods and the choice of filters to a periodic transonic flow 
with a shock. 

Proofs of the stability and convergence of Chebyshev approximations have 
the added complication of the boundary conditions and the weight, w(x), which 
is unbounded at the endpoints. In particular, the case where a(x) changes 
sign makes it difficult to show stability. Gottlieb [18] has proved stability 
of the straightforward Chebyshev collocation for the special cases where a(x) 
» ±x. 

To show stability of Chebyshev approximations in general, the skew- 
symmetric form of the equations is needed. We will survey the convergence 
theory of Canuto and Quarteroni [23] for the special case of the hyperbolic 
boundary value problem with b(-l) > 0 
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+ (bu)^ + b^u = f x€(-l,l), t€(0,T] 

u(-l,t) * 0 te(0,T] 

u(x,0) = Uq xe(-l,l). 

The further assumtion is added that 


(59) 


1/2 b^ + b^ - 1/2 bw^w"^ > 0 for x€(-l,l). (60) 

For the Chebyshev weight w(x) = (1 - x^)”^/2^ ^gg gf integration by 
parts to get an energy estimate will give an unbounded boundary term evaluated 
at X - +1 (see Gottlieb and Orszag [!])• This has led to the use of a modi- 
fied weight and norm with which to prove stability and convergence. Let the 
new weight be w*(x) * (1 - x)w(x) so that w*(l) * 0. Then the following 
convergence estimate holds: 


Theorea 17: Suppose that u€L"(0,T; H^*(-1,1)) and b € H®^(-1,1) 

w w 

for 6 > 2, Then the skew-symmetric Chebyshev approximation to (59) 
satisfies 


iu - 


u II 
c . 


'(L^*) ” 
w 


CN 


2 - 6 , 


w 


)• 


(61) 


Proof : See Canute and Quarteroni (23], Theorem 2.3. Note: Their 

theorem actually allows for more general boundary conditions than we have men- 


tioned here 
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For computations which are not done in skew-symmetric form, such conver- 
gence estimates are not available in the general case • As indicated above , 
Gottlieb [18] has shown stability in some particular cases* However, Reyna 
[24] has shown that if b(x) is not strictly positive in the interval that a 
straight-forward Chebyshev approximation need not be stable. To stabilize the 
solutions he proposes the use of filtering. It is not sufficient, however, to 
simply smooth the Chebyshev coefficients. Rather, he shows that stability can 
! be proved if Legendre coefficients are computed from the Chebyshev ones, the 

Legendre coefficients are smoothed, and then transformed back. 

The stability and convergence of Chebyshev approximations to the hyper- 
bolic initial-boundary-value problem for systems 

u = Au -1 < X < 1, t > 0 (62) 

t X — — _ 

where u is an n-vector and A is a constant matrix has recently been proved 
by Gottlieb, Lustman, and Tadmor [25], [26]. Because this system is hyper- 
bolic, the matrix A can be assumed to have been diagonalized to 

0 

A = 

.0 A 

where A^ < 0, A^^ > 0 are diagonal matrices. 

Boundary conditions for which this system is well posed are of the form 

u^(-l,t) = Lu^^(-l.t) + g^(t) 
u^^(l,t) = R u"(l,t) + g^^(t) 



( 63 ) 
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where and represent the partition of u into inflow and outflow 

components (see Kreiss and Oliger [5]). Under the assumption that the bound- 
ary conditions are dissipative, the standard Chebyshev collocation is stable: 

Theorem 16: Under the assumption that |r| • |l|^ 1-6<1, the 

Chebyshev collocation method is stable for the system (62) with boundary con- 
ditions (63) in the sense that there exists a weighting pair w(x) and con- 
stants q and Uq 0 such that for all s with Re s = n > 

(h - nQ)llu^(x,s)ll^ < CN^‘^|g(s)|^ 

where u^ and g are the Fourier transforms of u^ and g. 

Proof: See Gottlieb, Lustman, and Tadmor [25]. 


2. SOME APPLICATIONS OF SPECTRAL COLLOCATION METHODS 

In this section, some recent developments in the application of spectral 
methods to problems in fluid mechanics are surveyed. Much current emphasis 
has involved making spectral methods more efficient and more applicable to 
problems with complicated geometries. This has lead to the development of 
spectral multidomain methods which eliminate the need for global mappings and 
to the development of iterative techniques for the rapid inversion of the full 
matrices which occur when implicit time discretizations are used. 
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2.1. METHODS FOR ELLIPTIC AND PARABOLIC PROBLEMS 

2.1.1. SPECTRAL MULTID(X1AIN METHODS 

Spectral nrultidomain methods have been developed in order to avoid the 
need for global mappings required by spectral methods in problems with compli- 
cated geometries. A complicated domain can be subdivided into several subdo- 
mains and individual spectral discretizations can be applied to each subdo- 
main. For elliptic and parabolic problems, for handling the interfaces, early 
work considered explicit enforcement of continuity (e.g. Orszag [27] and 
Morchoisne [28]). More recently, spectral element discretizations and en- 
forcement of global flux balance have been used. The spectral element methods 
retain the accuracy of spectral methods in the context of a geometrically 
flexible finite element formulation. Global flux conservation has been used 
effectively when the mappings and/or domain sizes vary radically across inter- 
faces. 

Consider first the solution of the (second-order, self-adjoint, elliptic) 
Helmholtz equation, 

7^u - X^u = f in D (64) 

with Dirichlet boundary conditions on the domain boundary, 3D. Following 
the lead of finite element techniques, the spectral element algorithm [29, 30] 
proceeds by recognizing the equivalence of (64) to maximization of the follow- 
ing variational form. 


-7u»Vu/2 - X^u^/2 - uf}dx. 


max 

U6H^ 


(65) 
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The variational form, (65), is preferable over the differential statement, 
(64), in that it requires less continuity of candidate solutions. 

The spectral element discretization proceeds by breaking up the computa- 
tional domain, D, into general quadrilateral elements. Within a given element 
k, the solution, geometry, and data are then expanded as tensor product 
Lagrangian interpolants through a set of specified collocation points. For 
instance, in two space dimensions, the solution u in element k is 
represented as, 

u*^(r,s) = I I h. (r)h (s) (66a) 

i j ^ J 

h (z ) = 6 , (66b) 

m n mm 


where r and s are the local elemental coordinates, the h^ are the 
Lagrangian interpolants, the z„ are the collocation points, and 6 is 

the Kronecker delta symbol. All summations run from 0 to N, where N is the 
order of the Lagrangian interpolants in each element. 

The expansions (66a) are then inserted into (65), and the functional 
rendered stationary with respect to arbitrary variations in the nodal values, 
u^j , Direct stiffness summation [31] (which recognizes that the global 
approximation space must be C^) is then used to assemble the elemental equa- 
tions into the system matrix. It should be noted that, as regards the treat- 
ment of elliptic and parabolic equations, the "spectral element" recipe pre- 
sented here is very similar to earlier "p-type finite element" methods [32] 


and the "global element" method [33]. 


It is clear from the above representation, (66), that the global inter- 
polant space is only C^, that is, that the approximation space suffers dis- 
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continuities in derivative at elemental boundaries. Although this may appear 
to violate the basic smoothness required of spectral methods, this is not the 
case due to the fact that the variational formulation, (65), is used rather 
than the (unintegrated) Galerkin weighted-residual form. In particular, in 
the absence of "variational crimes", the spectral element method can be shown 
to achieve exponential convergence to smooth solutions as N, the order of 
(fixed) elements, is increased. For nonsmooth solutions (e.g., corner-induced 
singularities), high-order convergence is more difficult to obtain, however 
refinement techniques have been developed for the p-type finite element method 
[32], 

Variational crimes take the form of numerical quadrature errors and in- 
terpolation of boundary data. (Nonconforming elements are not considered.) 
In order to insure that these errors do not dominate the approximation errors, 
it is Important to correctly choose the collocation points of the Lagrangian 
interpolants. Earlier work on spectral elements used the Chebyshev colloca- 
tion points, as they are simple to evaluate and amenable to fast transform 
techniques. However, as the variational formulation (65) has essentially a 
unity weighting, it appears that a better choice is the Legendre-Lobatto 
points from the point of view of accuracy and efficiency of numerical 
quadratures [6, 34]. Although Legendre polynomials are less convenient than 
Chebyshev polynomials, are subject to round-off errors for high-order expan- 
sions, and cannot be "fast transformed", none of these objections are 
particularly oppressive for the relatively low-order expansions used in 
spectral element methods. 

As an example of the accuracy of a Legendre spectral element [34] (see 


[29, 30] for extensive discussion of Chebyshev-based techniques), consider the 
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problem 

V^u = 0 In D (67a) 

where D Is the domain defined in Figure 1, x€[0, 1], y€[0, 1 ^ 

Dirichlet boundary conditions are imposed such that the solution to the prob- 
lem is given by, 

u(x,y) = sin(x)e~y. (67b) 

The L error for the spectral element mesh shown in Figure 1 is plotted in 
00 

Figure 2, As expected from the analytic nature of the solution (67b) in the 
complex plane, exponential convergence is achieved as the order of the 
elements is increased* 

As another example of elliptic problems, consider the moving-boundary 
Stefan problem [34], given by 


V^e =0 in Dj, D 2 (68a) 

0 = 1/2 on (68b) 

70»n|_ + 3V0‘n|^ = 0 on 3D^ (68c) 

V0*n = 0 on 3 Dq (68d) 

0=1 + 1/2 cos2tix on 3D^ (68e) 

0 = 0 on 3D2» (68f) 


-33- 


where D^, D 2 > defined in Figure 3. Here the 

evaluation of _ and ^ refer to the Dj and D 2 sides of 9Dj, 
respectively. In point of fact, the time-dependent (parabolic) version of 
(68) was solved, approaching the steady-state only as t since the 

solution of parabolic equations involves at each time step the solution of an 
elliptic equation of the form (64), this aspect of the problem does not 
warrant separate discussion. 

Solution of the Stefan problem (68) illustrates several aspects of the 
spectral element method. First, since the interface is unknown and 

general, it demonstrates the ability to handle complex geometry. Second, 
though the solution suffers a discontinuity at the method has the 

ability to resolve certain non-homogeneities without losing "spectral 
accuracy". Figure 4 shows the interface position obtained with a Legendre 
spectral element method using a two-element mesh. In Figure 5, the associated 
temperature (9) distribution is given. High accuracy can be achieved with 
very few points. 

It is critical that the spectral element schemes not only be accurate, 
but also efficient as regards work required for a given level of accuracy , 
The key to the computational efficiency of the techniques is the sum factori- 
zation which follows from the tensor product representation, (66) , For in- 
stance, a typical elemental term in a two-dimensional Chebyshev spectral 
element equation is of the form, 



(69) 


where hj^, are defined as in (66), and all subscripts range from 0 to 
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N, the order of the polynomial space in each co-ordinate direction. Naive 
evaluation of this sum gives an operation count of O(N^), and O(N^) in 
three dimensions. This sum factorization is at the heart of both direct 
solvers using static condensation and fast eigenfunction solvers [35] and 
iterative solvers using conjugate gradient algorithms [36]. 

Another approach to handling domain interfaces was taken by Macaraeg and 
Streett [37], [38]. Within subdomains, the usual collocation procedure 

described in Section 1.2.2 is used. The interface values are computed by 
requiring that the solution be continuous and that the global flux be 
balanced. As an example of the procedure, consider the equation 

G(u) = F (u) - vu “ S(u) 

X XX 

u(-l) = a (70) 

u(l) = b 

where an interface is placed at x = x^* Integration of (70) from -1 to 
1 and the requirement that the jump in the flux [G] be zero at the interface 
yields 


*1 . 1 
/ S(u)dx = G(u)^j - / S(u)dx. 
-1 x^ 


(71) 


Numerical experiments show that spectral accuracy is retained. In two dimen- 
sions, the method has been used to solve Laplace's equation with discontinuous 
boundary conditions. 
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2.1.2 ITERATIVE SPECTRAL METHODS 

For evolution problems, explicit time-stepping can be extremely ineffi- 
cient. This is because the typical time-step limitation for spectral methods 

o / 

is proportional to l/N"^ for the advection equation and 1/N^ for the dif- 
fusion equation (where N is the number of modes) [39]* Hence, implicit time- 
stepping becomes a necessity. This results in a set of algebraic equations 
which are, in general, amenable to Iterative solution techniques only. Also, 
elliptic equations governing practical problems virtually require implicit 
iterative techniques. Since the condition number of the relevant matrices are 
large, preconditioned iterative schemes including multigrid procedures are the 
attractive choices. In this section, the fundamentals of iterative spectral 
methods are discussed with reference to an elementary example. 

For the purpose of illustration, consider the equation, 

= f, (72) 

periodic on [0,2 tt]. For the Fourier method, the standard choice of colloca- 
tion points is given in Equation (11). 

The Fourier collocation discretization of the equation (72) may be 
written 

LU = F, (73) 

where U = (u^, u^, F = (Fq* •••> ^ = C“^DC. 

Here C is the discrete Fourier transform operator, C ^ the inverse trans- 
form, and D the diagonal matrix denoting the first derivative operator in the 


Fourier space. Specifically, 
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and 



-2iTlk 

e 


(i-N/2) 

N 

) 



Kj - N/2) 
- 0 


j,k = 0, 1, N-1 

for J * 1,2, ..., N-1 
for j = 0 


(74) 


The eigenvalues of L are X(p) = ip, p = -N/2 + 1, ..., N/2 - 1, and the 
largest one grows as N/2. A preconditioned Richardson iterative procedure for 
solving Eq. (73) is 

V V + u)H“^ (F - LV) (75) 


where the preconditioning matrix, H, is a sparse, readily invertible approxi- 
mation to L. An obvious choice for H is a finite difference approximation 

to the first derivative. With the various possibilities for the 

eigenvalue spectrum of is given in Table I. Apparently, the staggered 

grid leads to the most effective treatment of the first derivative. This kind 
of preconditioning was successfully used in the semi-implicit time-stepping 
algorithm for the Navier-Stokes equations discussed in section 2.2 on Navier- 
Stokes Algorithms. The eigenvalue trends of that complicated set of vector 
equations are surprisingly well predicted by this extremely simple scalar 
periodic problem. 

Next, consider the second order equation 

-u = f on [0,2ir] (76) 

XX 


with periodic boundary conditions. A Fourier collocation discretization of 
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this equation is the same as Eq, (73) except for the diagonal matrix D which 
represents now the second derivative operator in the Fourier space. 


Djj = j - (j - f)^, j = 1, 2, N-1 


(77) 


(0, j = 0 

The eigenvalues of L are X(p) = p^ » P “ -N/2 + 1, ...» N/2-1. To make 
the case for the multigrid procedure (consisting of a fine-grid operator and a 
coarse— grid correction) as a preconditioner, assume H to be the identity 
matrix I in the iterative scheme (75). The iterative scheme is convergent if 
the eigenvalues, (1 - mX), of the iteration matrix [I - wL] satisfy 


1 - wX I < 1 . 


Each iteration damps the error component corresponding to X by a factor 
v(X) = 1 1-wX I . The optimal choice of X is that which balances damping of 
the lowest-frequency and the highest-frequency errors, i.e.. 


This yields 


(1 - o)X ) = - (1 - o)X ^ ) 
max min 


(I) 


SG 


TX + X . ) 
max min 


(78) 


(79) 


and the spectral radius 


^^max ^rain^ 

^SG Tx + X . ) * 

max min 


(80) 


9 2 

In the present instance, X ^^ ^ = N /4, X^^^ = 1, and thus = 1 - 8/N « 
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2 

This implies order N iterations are needed for convergence. This poor 
performance is due to balancing the damping of the lowest frequency eigenfunc- 
tion with the highest-frequency one. The raultigrid procedure exploits the 
fact that the lowest-f requency modes (|p| < N/4) can be damped efficiently on 
coarser grids, and settles for a relaxation parameter value which balances the 
damping of the mid-frequency mode (|p| = N/4) with the highest-frequency one 
(|p| = N/2). Table II provides the comparison of single-grid and raultigrid 
damping factors for N=64. The high frequencies from 16 to 32 are damped 
effectively in the multigrid procedure, whereas the frequencies lower than 16 
are hardly damped at all. But then some of these low frequencies (from 8 to 
16) can be efficiently damped on the coarser grid with N=32. Further coarser 
grids can be employed until relaxation becomes so cheap that all the remaining 
modes can be damped. In concrete terms, the ingredients of a multigrid tech- 
nique are a fine-grid operator, a relaxation scheme, a restriction operator 
which interpolates a function from the fine grid to the coarse grid, a coarse- 
grid operator, and a prolongation operator interpolating a function from the 
coarse grid to the fine grid. The fine grid problem for the present example 
may be written 

Lfyf « pf ^ (31) 

Let denote the fine-grid approximation. After the high-frequency content 

of the error has been sufficiently damped, attention shifts to the 

coarse grid. The coarse-grid problem is 

( 82 ) 


where 
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= R [F^ - V^] , 

R being the restriction operator. After a satisfactory approxmation, is 

obtained; the coarse-grid correction - RV^) is interpolated onto the 

fine grid by the prolongation operator P, yielding the corrected fine-grid 
solution 

V^-«- + P (V*^ - RV^) (83) 

I 

i 

The details of spectral multigrid techniques are furnished in [40]. Spectral 
multigrid techniques have been used to solve a variety of problems including 
the transonic full potential equation [41, 42 ]• Additional applications of 
spectral methods to compressible flows are described in [42]. 

[ 2.1.3 Convection-Dominated Flows 

A model for convection-dominated flows is the viscous Burger^ s equation, 

\ 9 

+ (u^) /2 = vu u(x,t=0) = -simrx, (84) 

t X XX ’ 

with boundary conditions u(-l) = u(l) =0, and "small dif f-usivlty," 

V = .01/ir [43]. The solution to this problem develops a near shock. This 

near shock is characterized by the time at which the derivative at the origin 
attains a maximum value, t^ax» value of its maximum derivative, 

|3u/9x| . The convective term is clearly dominant for short times, however 

the diffusion terra insures that the solution will be smooth. This convection- 
diffusion balance is a good model for the kind of phenomena that arise in 
solution of the Incompressible Navier-Stokes equations. The critical numer- 
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ical issues are numerical dispersion and diffusion. The former leads to in- 
correct propagation speeds of the shock affecting The latter leads to 
smearing of the shock affecting |3u/3x| . 

This problem has been solved by a variety of methods, including the spec- 
tral element method [43] and the explicit flux balancing method [37]. The 
spectral element calculations have used Crank-Nicolson in time on the diffu- 
sion terra and the resulting Helmholtz equation in space was solved using the 
variational methods presented in Section 2.1.1. The convective term was 
handled by explicit third-order Adams-Bashforth. Four elements were used 
covering the intervals [-1., -0.05], [-0.05, 0], [0., 0.05], [0.05, 1.] which 
cluster points around the location of high function variation. Macaraeg and 
Streett [37] used three subdomains with their flux conservation method. Table 
III presents a comparison of various methods on this model problem. 


2.2. INCOMPRESSIBLE NAVIER-STOKES EQUATIONS 

This section is devoted to a description of algorithms for the solution 
of the incompressible Navier-Stokes equations in primitive variable form. The 
algorithms are based on methods discussed in the previous section in the 
simplest context. For example, simulation of instability and transition to 
turbulence in a flat-plate boundary layer have used iterative methods 
described in section 2.1.2. The spectral element method has been used for a 
variety of flow computations, including the problem of flow past a cylinder. 
The Navier-Stokes equations in the so-called rotation form are 

q^ = qxoj+V • (y7q) - VP 


in D 



-Al- 


and 


V • q = 0 
q(x,0) = 


q = g 


in D (85) 

in D 

on 3D 


where q = (u,v,w) is the velocity vector, u = 7 x q the vorticity, 

P = p + 1/2 |q| the total pressure, y the variable viscosity, D the in- 

terior of the domain, and 3D its boundary. In the stability and transition 
problems, the domain D is cartesian and semi-infinite; periodic in the two 
horizontal directions (x,z), and bounded by a wall at y=0. Fourier colloca- 
tion can be used in the periodic directions (x,z) and Chebyshev collocation is 
used in the vertical (y) direction. The collocation points in the periodic 
directions are given by a relation similar to Eq. (11). The vertical extent 

of the domain 0 < y < » is mapped onto -1 <5 < +!• The velocities are 

defined and the momentum equations enforced at the points 

C = cos (^), j * 0, 1, ..., N (86) 

J y y 

The pressure is defined at the half points 


c J = cos[I-iiil^ ] , j = 0, 1, ..., N -1 (87) 

J+2 ^ 

and the continuity equation is enforced at these points* The staggered grid 
avoids artificial pressure boundary conditions , and precludes spurious pres- 


sure modes 
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After a Fourier transform in x and z, the temporal discretization 
(backward Euler for pressure, Crank— Nicolson for normal diffusion, and third 
or fourth-order Runge-Kutta for the remaining terms) of Eq. (85) leads to 

[I - MDM] Q + At Aq Vn = (88) 

- A^ 7 • Q = 0 

where 


,~n+l ~n+l ~n+l. 

Q “ {^Q » » • • • > } 


„ _ r~n+l ''n+1 ''n+1 

n - iPj/2 » P3/2 ‘ ***’ %-l/2^ 


V = {ik ik } 

X 3y^ X 


(89) 


M is the Chebyshev derivative operator, D the diagonal matrix with l/2yAt 
as its elements, and Ag is the interpolation operator from the half points 
to cell faces, vice versa. Obviously, the equations for each pair of 

horizontal wave number (k^ , k^) are indepenent and they can be written as 
the system 


LX = F 


(90) 


where X = [Q, n]. The iterative solution of this equation is carried out by 
preconditioning the system with a finite difference approximation on the 
Chebyshev grid, and applying a standard iterative technique such as 
Richardson, minimum residual or multigrid [44]. 
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The method described above solves the implicit equations together as a 
set* The extension of this method to the more general cases of interest such 
as those involving two or more inhomogeneous directions is not straightfor- 
ward* An alternative is the operator-splitting technique or the fractional 
step scheme [45]* This method yields implicit matrices which are positive 
definite and are easily amenable to iterative methods* In the first step, one 
solves the advection-diffusion equation 


^ ^ ^ A 

q =q xo) + V» (yVq ) 


(91) 


subject to the initial and boundary conditions 


q (x,t^) = q(x,t^). 


* * 

q = g on 3D* 


(92) 


Note that g has yet to be defined. In the second step, one solves for the 
pressure correction 

(93) 

(94) 

subject to the conditions 


** ** 

** 

V • q =0 


** * * * . 
q (x,t ) = q (x,t ) 




n = g • n 


in D 
in 


(95) 


where n is the unit normal to the boundary 


Further, the tangential 
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component of the Eq. (95) holds on the boundary, i.e., 

. T = - VP . X In 3D (96) 

A ?Hr 

where x is a unit tangent vector to the boundary. Now g is defined 
[45] as (using Taylor expansion in t) 

g* • n = (g” + At g”) • n (97) 

g • X = [g*' + At (g” + VP”)1 • X . 

Eq. (91) is discretized in the usual spectral collocation manner. After a 
temporal and spatial discretization of Eq. (93), the boundary conditions are 
built into the relevant matrix operators, and then a discrete divergence is 
taken. This results in a discrete Poisson equation (with aS many algebraic 
equations as unknowns) for pressure, which can be solved by standard iterative 
techniques including the multigrid method. 

Spectral element methods have been applied to the incompressible Navier- 
Stokes equations (85). In addition to (98), the uncoupled (passive) or 
coupled (natural convection) energy equation is also often of interest. The 
time discretization used for the Navier-Stokes equations is either a Green's 
function technique [29] or an operator splitting scheme [30]. Both of these 
methods reduce (85) at each time step to an initial convective step, followed 
by a Stokes problem consisting of a sequence of Poisson and Helmholtz equa- 
tions. The spatial discretizations discussed above in Section 2.1.1 are 


directly applicable to these Navier-Stokes subproblems 
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Spectral element methods have been applied to the simulation of numerous 
flows in the Reynolds number range 0 < R < 1500 [36, 46 - 47]. An example is 
provided by the classical problem of flow past a cylinder. Results are pre- 
sented here for R = 100, based on freestream velocity and cylinder diameter, 
for times sufficiently large that the flow has reached a steady-periodic 
state. Figure 6 shows the spectral element mesh used, and Figures 7 and 8 
show the streamlines and isotherms, respectively, at one time in the periodic 

flow cycle. The thermal boundary conditions are T = T far from the 

00 

cylinder, and T = on the cylinder surface. The isotherm pattern clearly 
reveals the spatial structure of the von Karraan vortex street. Note the mini- 
mal numerical dispersion in the scheme, as evidenced by the clear identity of 
the shed packets of fluid and the absence of trailing waves in the wake. More 
details of these cylinder calculations, as well as comparisons with previous 
numerical work and experiment, can be found in [36]. 


2 . 3 HYPERBOLIC EQUATIONS 

Here, the application of spectral methods to the solution of invlscid 
compressible flow problems is surveyed. Methods for such problems are not 
nearly so advanced as those for incompressible flows. The survey is limited 
to methods for the solution of the Euler equations of gas-dynamics governing 
some flows of aerodynamic interest. For the solution of the full potential 
equation for transonic flows, see Streett, et al. [42]. 

The Euler equations of gas-dynamics are a coupled system of nonlinear 
hyperbolic equations which (in one dimension) are usually written in the con- 


servative form 
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ut + Fx(u) = H. (98) 

Typically, spectral discretization in space and explicit finite difference 
discretization in time are used. The discontinuous solutions of this set of 
equations have been obtained in the case of a shock tube (Gottlieb, Lustman, 
and Orszag [48], Cornille [49]), quasi-one-dimensional flow in a nozzle (Zang 
and Hussaini [50]) and for the astrophysical problem of shocked flow in a 
galaxy (Hussaini, et al, [51]), 

The astrophysical problem is the most challenging one-dimensional com- 
pressible flow problem for which shock capturing has been attempted with a 
Fourier spectral method. It contains a strong shock and an adjacent strong 
expansion. Unlike problems with weak shocks and expansions, it was necessary 
to apply strong filtering to stabilize the numerical solution. The result of 
this drastic filtering was a reduction of the order of accuracy# Even in the 
smooth parts of the solution away from the shock, the accuracy was only first 
order. In view of the extra work involved to compute the spectral approxima- 
tions, it is not clear that spectral methods with filtering are a viable 
alternative to finite difference methods when strong shocks are captured. 

An alternative to capturing shocks is to treat them as boundaries. In 
this case, it is possible to compute the solutions using the nonconservative 
form 

\ = E (99) 

along with an ordinary differential equation for the motion of the shock, A 
number of two dimensional shock-fitted solutions are described in Hussaini, et 
al. [52], These solutions include a shock/ turbulence interaction. 
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Shock/ vortex interaction and supersonic flow past a cylinder* When shocks are 
fitted, spectral methods do indeed outperform typical second order finite 
difference methods, as long as the solution is adequately resolved. Kopriva, 
et al. [53] compared the performance between MacCormack's method and the spec- 
tral collocation method for the shock/plane wave interaction problem and for 
the Ringleb problem. A comparison of the accuracy of the finite difference 
method vs. the spectral method is shown in Table IV. 

A multidomain method for the nonconservative form of the Euler equations 
suitable for use with shock-fitting has been described by Kopriva [54]. In 
each subdomain, the usual collocation method (Section 1.2.2) is applied. At 
interfaces, however, a weighted average of the derivatives is used. In one 
dimension, 

= E (100) 

where denotes the solution vector at the interface and the derivatives 

superscripted with the L and R denote the left and right computed spectral 
approximations. For consistency, A^ + A^ = A. The weighting corresponds to 
an upwind approximation 

A^ = 1/2(A + |a|) = 1/2(A - |a|) 

where |a| = zIaIz"^ and Z is the matrix of right eigenvectors. For many 

applications, this can be simplified by replacing |a| by a diagonal approxi- 

mation IaI » X I where |X i . < X < Ixl is an approximation to the 

' ' ' 'min ' 'max 

o *f cron-ir o 1 HOC nf A 

W A.A. ▼ AX# 
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Subdividing the domain retains spectral accuracy. Table V shows the 
performance of a two domain computation of a one dimensional 2x2 
hyperbolic system 

= 0 (102) 

X 

from [54] with solution u and v for an equal number of points on each side 
of the interface. 

For a given number of grid points, it is possible to obtain solutions 
with a multidomain spectral method which are significantly better than the 
single domain method. For the two dimensional nonlinear Ringleb problem 
computed by Kopriva [54], Table VI shows the effect on the error of a four 
domain division in which the position of the streamwise interface is varied. 
The accuracy is best when rapid changes in the solution are best resolved. 

The shock-vortex interaction problem described by Kopriva [55] provides 
another example of the advantage of a multidomain method over a single domain 
method. A two dimensional region between the shock and an arbitrary upstream 
boundary is mapped onto a square. The shock moves downstream where it en- 
counters a vortex. The interaction of the shock and the vortex creates a 
circular sound wave centered on the vortex. Because the physical domain is 
continually increasing in size, the resolution of the solution decreases with 
time. 

The single domain solution to the shock-vortex problem cannot be computed 
without added smoothing. Figure 9a shows the pressure contours with no 
smoothing. The numerical oscillations in the pressure are of the same order 
as the sound pressure wave created by the interaction. If the region between 


u 


1 2' 


u' 


+ 




V .. 


.2 1. 


V 
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the shock and the upstream boundary Is subdivided and some of the subdomains 
are allowed to move with the shock, smoothing is not required. Figure 9b 
shows the pressure contours of a two domain calculation with the same number 
of grid points in the horizontal direction. The horizontal numerical oscilla- 
tions are no longer present and the sound pressure wave is clearly visible. 
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Figure 1 . 


Figure 2. 


Figure 3. 


Figure 4, 


Figure 5. 


Figure 6. 


Figure 7 


Figure 8 


Figure 9 


FIGURE CAPTIONS 


The domain and Legendre spectral element discretization used for 
solution of the Laplace equation described in the text. Although 
results are given for a single isoparametric element, similar re- 
sults are obtained with multiple elements. (Legendre results due 
to E. M. Ronquist.) 


A plot of the L error as a function of the number of points 
in one direction for solution of Laplace's equation in the domain 
shown in Figure 1 . 


Description of the domain and boundaries for the Stefan problem 
presented in the text. (Stefan problem results due to E. M. 
Ronquist.) 


Spectral element prediction for the position of the interface, 

9Dj, for the Stefan problem described in the text. The spectral 
element mesh uses two elements, one each in each phase. 


Temperature (0 ) distribution for the Stefan problem described 
in the text. Note the discontinuity of slope at the interface, 
9Dj. 


Spectral element mesh used for simulation of flow past a cylin- 
der. Note the flexible resolution afforded by the elemental 
decomposition. (Flow past a cylinder results due to G. E. 
Karniadakis.) 


Instantaneous streamlines of the cylinder flow at a Reynolds 
number of R = 100. 


Instantaneous isotherms of the cylinder flow at a Reynolds number 
of R = 100 (Prandtl number of unity). The von Karman vortex 
street can be clearly seen in the temperature distribution behind 
the cylinder. 


Pressure contours for a shock/ vortex interaction. (a) Single 
domain calculation, (b) Three vertical domain calculation. 
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Figure 7. 



Figure 8 



64 



Figure 9a 




Figure 9b 


Table 1. Preconditioned Eigenvalues for One-dlaenslonal First Derivative 

Model Problem 


Preconditioning 


Eigenvalues 


Central Differences 


kAx 

sln(kAx) 


High Mode Cutoff 


kAx 

sin(kAx) 


IkAxI < (2x/3) 


(2x/3) < IkAxI < x 


One-sided Differences 


-l(kAx/2) kAx/2 
' sin((kAx)/2) 


Staggered Grid 


(kAx) /2 
sin( (kAx)/ 
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Table Ills Goiqiarlson of Hethods for Solution of 
Burger' a Equation (frou Ref. 43, 37) 


Method 

Interval 


TT» t 

max 

N/M 

At»Tf 

- Fourier Galerkln 

t-1,1] 

151.942 

1.6035 

682/1024 

5.10-^ 



142.665 

1.60 

682/1024 

1 

o 



148.975 

1.603 

170/256 

1 

o 

• 



142.313 

1.60 

170/256 

10“2 

- Fourier Pseudospectral 

1-1,1] 

142.606 

1.60 

256/256 

10"2 



144.237 

1.60 

128/128 

10”2 

- ABCN collocation 

[-1.1] 

145.877 

1.60 

512 

5.10"^ 

+coordinate transform 

[-1,1] 

152.123 

1.60 

64 

10~2 

- Spectral Element 

i 

[-1.1] 

152.04 

1.6033 

16 X 4 

10"2/6 

- FD 

[-1,1] 

150.1 

1.63 

81 

10-2 

- Chebyshev 

I 





ABCN spectral 


152.05 

1.60 

64 

1/300 

Rosenbrook spectral | 


151.998 

1.60 

64 

10“2 


[0,1] 

150.144 

1.60 

32 

10"2 

ABCN collocation 

[0,1] 

152.126 

1.60 

64 

10“2 

- Flux balance 

[-1,1] 

152.00011 





Analytical 


152.00516 1.6037 
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Table IV. Maxiam &ror la p for HacCoraack aod Spectral Coaputatioa 

of Transonic Rlngleb flow 


Grid 

MacCormack 

Spectral 

9x5 

CM 

1 

o 

X 

sO 

• 

CM 

2.2 X 10"2 

17 X 9 

1.1 X 10“2 

1.9 X 10"^ 

33 X 17 

3.2 X 10"3 

5.0 X 10”5 

Table V. Solutions 

to (102) with Equal Number of Points on Each Side 
of the Interface 

N 

Error in u 

Error in v 

8~ 

1.57 X 10"2 

1.49 X 10”2 

16 

4.15 X 10"^ 

4.86 X 10"^ 

32 

1.91 X 10"^ 

1.91 X 10“^ 


Table VI. Effect of Streaawlse Mesh Distribution on Slngleb Calculation 


Grid 


Division 


Maximum Error 


8+8 
8+8 
16 (SD) 
10 + 6 


0.45 + 0.55 
0.50 + 0.50 


7.8 X 10”^ 
9.3 X 10”^ 

1.9 X 10"3 


1.2 X 10 


-2 


0.34 + 0.66 
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